Joint AVO inversion, wavelet estimation and noise-level estimation using a spatially coupled hierarchical Bayesian model

The main objective of the AVO inversion is to obtain posterior distributions for P-wave velocity, S-wave velocity and density from specified prior distributions, seismic data and well-log data. The inversion problem also involves estimation of a seismic wavelet and the seismic-noise level. The noise model is represented by a zero mean Gaussian distribution specified by a covariance matrix. A method for joint AVO inversion, wavelet estimation and estimation of the noise level is developed in a Bayesian framework. The stochastic model includes uncertainty of both the elastic parameters, the wavelet, and the seismic and well-log data. The posterior distribution is explored by Markov-chain Monte-Carlo simulation using the Gibbs’ sampler algorithm. The inversion algorithm has been tested on a seismic line from the Heidrun Field with two wells located on the line. The use of a coloured seismic-noise model resulted in about 10% lower uncertainties for the P-wave velocity, S-wave velocity and density compared with a white-noise model. The uncertainty of the estimated wavelet is low. In the Heidrun example, the effect of including uncertainty of the wavelet and the noise level was marginal with respect to the AVO inversion results. I N T R O D U C T I O N Geophysical measurements are of crucial importance in making models of the structures and the physical properties of the subsurface. The geophysical measurements are often strongly affected by noise and measurement uncertainty, and the established subsurface models may be highly uncertain. Quantification of the uncertainty is important to appreciate these subsurface models correctly. In a Bayesian setting, available prior knowledge is combined with the information contained in the measured data (Tarantola and Valette 1982; Duijndam 1988a,b; Malinverno 2000; Scales and Tenorio 2001; Ulrych, Sacchi and Woodbury 2001). The prior knowledge about the model parameters is specified by probability density functions where the prior belief and the corresponding uncertainty are defined. The relationship between the model parameters and the measured data is described by the likelihood model. The solution of a Bayesian inverse problem is represented ∗E-mail: [email protected] by the posterior distribution, which provides both the most probable solution and information about the corresponding uncertainty. Amplitude versus offset (AVO) inversion can be used to extract information about the elastic subsurface parameters utilizing the angle dependency of the reflection coefficient (Smith and Gidlow 1987; Hampson and Russell 1990; Buland et al. 1996; Gouveia and Scales 1998; Wang 1999). The AVO inversion problem can be linearized if an appropriate processing sequence is applied to the seismic data prior to the inversion. Important elements in such processing are the removal of the moveout, multiples, and the effect of geometrical spreading and absorption. The seismic data should be prestack migrated (Wang, White and Pratt 2000; Buland and Landrø 2001) and transformed from offset to reflection angle. A Bayesian linearized AVO inversion is defined in Buland and Omre (2003a), where an explicit analytical form of the posterior distribution is derived under Gaussian model assumptions. The explicit analytical form of the posterior distribution provides a computationally fast inversion method suitable for inversion of large C © 2003 European Association of Geoscientists & Engineers 531 532 A. Buland and H. Omre 3D seismic data sets. However, the elastic parameters are not laterally coupled, so the inversion is performed independently for each bin gather. Furthermore, the seismic wavelet and the noise covariance are assumed to be known prior to the inversion. A Bayesian method for estimation of the wavelet and the noise covariance is presented in Buland and Omre (2003b). In this paper, a more realistic and complex statistical model is defined for the linearized AVO inversion problem. Firstly, a spatial coupling of the model parameters is imposed by a spatial correlation function. This ensures a spatially consistent solution. Secondly, the solution is obtained both from seismic prestack data and well logs. The spatial coupling of the model parameters is required for the integration of seismic and well-log data. Thirdly, the AVO inversion, the wavelet estimation and the estimation of the noise level are carried out simultaneously. The uncertainty of the wavelet and the noise level, and the corresponding effect on the estimated elastic parameters, are integrated in the algorithm. A general hierarchical Bayesian model is defined. This means that the statistical parameters specifying probability density functions are also regarded as uncertain, for example the mean vector and the covariance matrix for a Gaussian variable. These statistical parameters are assigned prior distributions specified by hyperparameters. For trivial Bayesian problems, analytical expressions for the posterior distributions can often be found. In the current case, no analytical solution exists, but the posterior distribution can be explored by Markov-chain Monte-Carlo (MCMC) simulation (Gilks, Richardson and Spiegelhalter 1996; Chen, Shao and Ibrahim 2000). The methodology is presented in the following sections, firstly in rather general terms, and then illustrated by an inversion example of a real data set from the Heidrun Field. M E T H O D O L O G Y An isotropic, elastic medium is described completely by the material parameters {α(x, t), β(x, t), ρ(x, t)}, where α, β and ρ are P-wave velocity, S-wave velocity and density, x is the lateral location, and t is the two-way vertical seismic traveltime. A weak contrast reflectivity function for PP reflections is (Aki and Richards 1980; Buland and Omre 2003a) c(x, t, θ ) = aα(x, t, θ ) ∂ ∂t ln α(x, t) + aβ (x, t, θ ) ∂ ∂t ln β(x, t) + aρ(x, t, θ ) ∂ ∂t ln ρ(x, t), (1) where θ is the angle of reflection, aα = (1 + tan2 θ )/2, aβ = −4(β/α)2 sin2 θ and aρ = (1 − 4(β/α)2 sin2 θ )/2. The inversion algorithm requires that aα, aβ and aρ are defined from a prior known background model. Motivated by the form of the reflectivity function in expression (1), let the unknown model parameter vector be m(x, t) = [ln α(x, t), ln β(x, t), ln ρ(x, t)], (2) where superscript T denotes transpose. The seismic data are represented by the convolutional model dobs(x, t, θ ) = ∫ s(τ, θ ) c(x, t − τ, θ ) dτ + ed(x, t, θ ), (3) where s is the wavelet, and ed is an error term. The wavelet is allowed to be angle dependent, but should be independent of the lateral location x. The wavelet is assumed to be stationary within a limited target window. The seismic data dobs(x, t, θ ) are available in a discrete form, denoted dobs. In general, the discretization of the subsurface parameters m(x, t) should be determined by the intrinsic variability and future use of the inverted parameters. An identical lateral and temporal sampling of the model parameters and the seismic data is often chosen. Usually, this is an adequate choice, but the methodology is not restricted to equal sampling. Let a discrete representation of m(x, t) be written m. Further, let s be a discrete representation of the seismic wavelet and let wobs be a discrete representation of the well-log data. The stochastic dependency model The variables and the problem structure are graphically displayed by a directed acyclic graph (DAG) in Fig. 1 (see e.g. Figure 1 The stochastic model represented by a directed acyclic graph. The nodes represent stochastic variables, the thin arrows represent probability dependencies, and the thick arrows represent deterministic relationships. C © 2003 European Association of Geoscientists & Engineers, Geophysical Prospecting, 51, 531–550 Joint AVO inversion, wavelet estimation and noise-level estimation 533 Lauritzen 1996; Spiegelhalter et al. 1996). All nodes in the graph represent stochastic variables, where μ and Σ denote expectation vectors and covariance matrices. The arrows represent the causal dependency structure, where thin lines indicate probability dependencies and thick lines represent deterministic relationships. The DAG can be interpreted as a family tree where the parent nodes point to its children. When probabilistic dependencies are considered, nodes connected by a deterministic link merge into a single node. The expectation vectors μw and μd are deterministically determined from m and from m and s, respectively. The parents of the well-log data are therefore Pa(wobs) = {m, Σw}, and the parents of the seismic data are Pa(dobs) = {m, s, Σd}, where the notation Pa(v) denotes the parents of a node v. Before any data are observed, the marginal distributions of the parents are independent, but when their child is observed, the properties of the parents become dependent. An example is m and s which are a priori independent, but become dependent when the observed seismic data dobs are included in the model. The likelihood model The error term in the convolutional model in expression (3) is assumed to be zero mean Gaussian with covariance Σd. The likelihood model for the seismic data dobs is then Gaussian, compactly denoted dobs |μd,Σd ∼ Nnd (μd,Σd), (4) where nd is the dimension of the seismic data vector dobs. A definition of the multi-Gaussian distribution is given in Appendix A. The expectation vector μd is deterministically determined by the convolution of the reflection coefficients computed from m with the wavelet s. The relationship between the well observations wobs and the total model parameter vector m can be written wobs = Pm + ew, (5) where P is a design matrix of zeros and ones which defines the locations of the wells relative to m, and ew is an error term. If we assume that the well-log error ew is zero mean Gaussian with covariance Σw, then the likelihood model for the well-log information is Gaussian, wobs |μw,Σw ∼ Nnw (μw,Σw), (6) where nw is the dimension of wobs and μw = Pm. The well locations relative to the seismic data are assumed to be correct. Also the conversion of the well logs from depth to seismic traveltime is assumed to be correct. In principle, this uncertainty related to the depth-to-time transform of the well logs could have been included in the present method by simulation of possible shift, stretch and squeeze of the log time axis (Buland and Omre 2003b). The prior model The model parameters m are a priori assumed to be Gaussian, m |μm,Σm ∼ Nnm (μm,Σm), (7) where nm is the dimension of m, and μm and Σm are the expectation vector and covariance matrix, respectively. The expectation is usually constant or slowly varying for each of the parameters ln α, ln β and ln ρ in m. The covariance defines the variances and the correlations between the different elements in m. A spatial coupling of the parameters is imposed through the covariance matrix by a spatial correlation function. Also the wavelet s is a priori modelled with a Gaussian distribution, s |μs,Σs ∼ Nns (μs,Σs), (8) where ns is the dimension of s, and μs and Σs are the expectation vector and the covariance matrix, respectively. The length of the wavelet is fixed. In general, the length of the wavelet should be as short as possible, but still long enough to represent the link between the model parameters and the seismic data as defined by the convolutional model. Expected smoothness of the wavelet can be imposed through the covariance matrix by a temporal correlation function. The DAG in Fig. 1 represents a hierarchical Bayesian model, where the expectations and covariances are stochastic. The prior expressions for m and s in expressions (7) and (8) can be regarded as the first level of a hierarchical prior model, the structural portion, while the prior distributions for the expectation vectors and covariance matrices form the second level, the subjective portion of the prior (see Robert 1994; Carlin 1996). A mathematically convenient class of prior distributions for the second level is the conjugate prior distributions (Robert 1994). These distributions have the same parametric prior and posterior forms, but with different parameters. The conjugate distributions for expectation and variance in a Gaussian model are the Gaussian and the inverse gamma distributions, respectively. The uncertainties of μm and μs are therefore modelled by Gaussian distributions, μm ∼ Nnm ( μμm ,Σμm ) , (9) C © 2003 European Association of Geoscientists & Engineers, Geophysical Prospecting, 51, 531–550 534 A. Buland and H. Omre μs ∼ Nns ( μμs ,Σμs ) , (10) where the expectations μμm and μμs , and the covariances Σμm and Σμs , are fixed hyperparameters. The expectations μw and μd are deterministically determined and do not need prior distributions. If the covariance matrices Σm, Σs, Σw and Σd are known up to unknown multiplicative variance factors, the covariances can be written Σi = σ 2 i Σ0,i , (11) where i is one of {m, s, w, d}. The uncertainty of the unknown variance factor σ i is modelled by the inverse gamma distribution, σ 2 i ∼ IG(γi , λi ), (12) where the hyperparameters γ i and λi define the prior expectation and variance (see Appendix B). If a covariance matrix is completely unknown, a conjugate prior distribution can be obtained using the inverted Wishart distribution. Then the conditional distribution given the corresponding sample covariance matrix is also inverted Wishart distributed (see Appendix C). The posterior model The objective of the inversion is to obtain the posterior distributions for the variables involved, based on the observed seismic data dobs and the well logs wobs. Let the vector v represent the unknown quantities, that is, all variables in the DAG in Fig. 1 except dobs and wobs. Further, let o be a vector containing the observed data, o = [ wobs dobs ] . (13) The posterior distribution for v given dobs and wobs can be written p(v | o) = p(v, o) p(o) , (14) where p(v, o) is the complete joint distribution. The pdf p(o) does not contain unknown variables, and can therefore be regarded as a constant. The complete joint distribution for a model represented by a DAG can generally be written (Spiegelhalter et al. 1996) p(v, o) = ∏ xi ∈{v,o} p(xi | Pa(xi )). (15) An expression for the posterior distribution is now defined from expressions (14) and (15), and the DAG in Fig. 1, p(m, s,μm,μs,Σm,Σs,Σw,Σd | dobs, wobs) ∝ p(dobs | m, s,Σd)p(wobs | m,Σw)p(m |μm,Σm) p(s |μs,Σs)p(μm)p(μs)p(Σm)p(Σs)p(Σw)p(Σd), (16) where the proportionality is caused by the unknown constant probability density functions involving dobs and wobs. The Gibbs’ sampler algorithm The posterior distribution can be explored by MCMC simulation. The Gibbs’ sampler is perhaps the best known and most popular of the MCMC algorithms. The name of the algorithm was introduced by Geman and Geman (1984) who worked with Gibbs’ distributions on lattices, but the name is misleading as the application of the algorithm is general, and not restricted to Gibbs’ distributions. A general presentation of MCMC can be found in Gilks et al. (1996) and Chen et al. (2000). The Gibbs’ sampler algorithm works on the complete joint distribution and requires full conditional distributions for each single variable given all the others. One iteration consists of drawing new samples for the unknown quantities. A new sample is drawn conditioned on the current state of the other elements in v. Once an element is drawn, it goes into the current state of v. The pseudo-code of the Gibbs’ sampling algorithm can be written Initiate: Set arbitrary v(0) where p(v(0) | o) > 0. Iterate: For i = 1, 2, . . . , draw v 1 ∼ p(v1 | v−1, o), .. v j ∼ p(vj | v− j , o), .. v(i) n ∼ p(vnv | v−n, o), where v− j = {v 1 , . . . , v j−1, v(i−1) j+1 , . . . , v(i−1) nv }. Instead of drawing single elements from v as shown above, a slightly more general algorithm is obtained by allowing for simultaneous drawing of a group of elements. An example is to draw the complete wavelet s instead of updating single samples of the wavelet. The elements in the model parameter vector m can also be grouped, for example into the complete vector m, or some partition of m. The full conditional distributions needed in the Gibbs’ sampling are specified in Appendix D. C © 2003 European Association of Geoscientists & Engineers, Geophysical Prospecting, 51, 531–550 Joint AVO inversion, wavelet estimation and noise-level estimation 535 I N V E R S I O N E X A M P L E O F H E I D R U N D ATA The Heidrun Field is an oil and gas field, offshore midNorway. Geologically, the field comprises a heavily faulted −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 A B 11 degrees